package scu.maqiang.tiangong;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.fes.FES3T42;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

import java.io.File;
import java.io.FileNotFoundException;
import java.util.Arrays;
import java.util.Scanner;
import java.util.function.Consumer;

public class CTE_2C {
    public static void main(String[] args) {
        //String dir = "C:\\Users\\Administrator\\OneDrive\\材料-马强\\DSMC表面热流压力计算数据\\";
        //String[] TG2CData = {dir + "TG-2C_80000m.dat"};
        //importData(TG2CData, null);
        computeCTE_2C();
    }

    public static double[][][] importData(String[] dataFiles, Mesh3T4 objectMesh) {
        int fileNumber = dataFiles.length;
        double[][][] allDatas = new double[fileNumber][][];
        String firstFile = dataFiles[0];
        Scanner scanner = null;
        try {
            scanner = new Scanner(new File(firstFile));
        } catch (FileNotFoundException e) {
            e.printStackTrace();
        }
        String title = scanner.nextLine();
        int equalIdx = title.indexOf("=");
        String variables = title.substring(equalIdx + 1);
        System.out.println(variables);
        String[] variable = variables.split(",");
        int numberOfVariables = variables.length() - 3;
//        System.out.println(variable.length + Arrays.toString(variable));
        String info = scanner.nextLine();
        int NequalIdx = info.indexOf("N=");
        int EequalIdx = info.indexOf("E=");
        String subString1 = info.substring(NequalIdx + 2);
        int DouIdx1 = subString1.indexOf(",");
        String strNumberOfNodes = subString1.substring(0, DouIdx1).trim();
        String subString2 = info.substring(EequalIdx + 2);
        int DouIdx2 = subString1.indexOf(",");
        String strNumberOfElements = subString2.substring(0, DouIdx2).trim();
        //System.out.println(strNumberOfNodes);
        //System.out.println(strNumberOfElements);
        int numberOfNodes = Integer.parseInt(strNumberOfNodes);
        int numberOfElements = Integer.parseInt(strNumberOfElements);
//        System.out.println(numberOfNodes);
//        System.out.println(numberOfElements);
        scanner.close();

        Scanner[] filesScanner = new Scanner[fileNumber];
        try {
            for(int i = 0; i < fileNumber; i++) {
                filesScanner[i] = new Scanner(new File(dataFiles[i]));
                filesScanner[i].nextLine();
                filesScanner[i].nextLine();
            }
        } catch (FileNotFoundException e) {
            e.printStackTrace();
        }

        int numberOfNodesInMesh = objectMesh.getNv();
        for(int i = 0; i < fileNumber; i++) {
            allDatas[i] = new double[numberOfVariables][numberOfNodesInMesh];
        }
        double[] xyz = new double[3];
        for(int i = 0; i < numberOfNodes; i++) {
            xyz[0] = filesScanner[0].nextDouble();
            xyz[1] = filesScanner[0].nextDouble();
            xyz[2] = filesScanner[0].nextDouble();
            for(int k = 1; k < fileNumber; k++) {
                filesScanner[k].nextDouble();
                filesScanner[k].nextDouble();
                filesScanner[k].nextDouble();
            }
            double min = 1.0e10;
            int idx = -1;
            for(int k = 0; k < numberOfNodesInMesh; k++) {
                double[] node = objectMesh.getNode(i);
                double error = MVO.distance(xyz, node);
                if (min > error) {
                    min = error;
                    idx = k;
                }
                if(error < 1.0e-4) {
                    for(int kk = 0; kk < fileNumber; kk++) {
                        for(int mm = 0; mm < numberOfVariables; mm++) {
                            allDatas[kk][mm][k] = filesScanner[kk].nextDouble();
                        }
                    }
                    break;
                }
            }
        }
        return allDatas;
    }

    public static void computeCTE_2C() {
        String dir = "C:\\Users\\Administrator\\OneDrive\\TianGongComputation\\TianGong2CComputations\\";
        String[] heatFluxData = {dir + "HeatFluxData_110000m.txt", dir + "HeatFluxData_107500m.txt", dir + "HeatFluxData_105000m.txt", dir + "HeatFluxData_102500m.txt",
                dir + "HeatFluxData_100000m.txt", dir + "HeatFluxData_97500m.txt",  dir + "HeatFluxData_95000m.txt",  dir + "HeatFluxData_92500m.txt",
                dir + "HeatFluxData_90000m.txt",  dir + "HeatFluxData_87500m.txt",  dir + "HeatFluxData_85000m.txt",  dir + "HeatFluxData_82500m.txt",
                dir + "HeatFluxData_80000m.txt"};
        String[] pressureData = {dir + "pressureData_110000m.txt", dir + "pressureData_107500m.txt", dir + "pressureData_105000m.txt", dir + "pressureData_102500m.txt",
                dir + "pressureData_100000m.txt", dir + "pressureData_97500m.txt",  dir + "pressureData_95000m.txt",  dir + "pressureData_92500m.txt",
                dir + "pressureData_90000m.txt",  dir + "pressureData_87500m.txt",  dir + "pressureData_85000m.txt",  dir + "pressureData_82500m.txt",
                dir + "pressureData_80000m.txt"};
        Mesh3T4 mesh2C = new Mesh3T4().readMesh(dir + "Whole2CMesh2.mesh", MeshFormat.mesh);
        mesh2C.toTecplot(dir + "My2CMesh.dat");
        int dataCount = heatFluxData.length;
        int nodesInMesh = mesh2C.getNv();
        double[][] hfDataInMesh = new double[dataCount][];
        double[][] prDataInMesh = new double[dataCount][];
        for(int i = 0; i < dataCount; i++) {
            hfDataInMesh[i] = MVO.readArray(heatFluxData[i]);
            prDataInMesh[i] = MVO.readArray(pressureData[i]);
        }
        double time = 0.0;
        double dt = 1.0;
        double[] timePoint = {0.0, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500};
        double endTime = timePoint[timePoint.length - 1];
        double[] previousFlux = new double[nodesInMesh];
        double[] previousPres = new double[nodesInMesh];
        double[] currentFlux = new double[nodesInMesh];
        double[] currentPres = new double[nodesInMesh];
        Consumer<Double> setCurrentFlux = currenttime -> {
            if(Math.abs(currenttime) < 1.0e-6) {
                Arrays.setAll(currentFlux, i -> hfDataInMesh[0][i]);
                return;
            }
            int location = 0;
            for(int i = 0, n = timePoint.length; i < n; i++) {
                if (currenttime > timePoint[i]) {
                    location++;
                }
            }
            int currentLocation = location;
            Arrays.setAll(currentFlux, i -> {
                double timep = (currenttime - timePoint[currentLocation - 1]) / (timePoint[currentLocation] - timePoint[currentLocation - 1]);
                return (hfDataInMesh[currentLocation][i] - hfDataInMesh[currentLocation - 1][i]) * timep;
            });
        };
        Consumer<Double> setCurrentPressure = currenttime -> {
            if(Math.abs(currenttime) < 1.0e-6) {
                Arrays.setAll(currentPres, i -> prDataInMesh[0][i]);
                return;
            }
            int location = 0;
            for(int i = 0, n = timePoint.length; i < n; i++) {
                if (currenttime > timePoint[i]) {
                    location++;
                }
            }
            int currentLocation = location;
            Arrays.setAll(currentPres, i -> {
                double timep = (currenttime - timePoint[currentLocation - 1]) / (timePoint[currentLocation] - timePoint[currentLocation - 1]);
                return (prDataInMesh[currentLocation][i] - prDataInMesh[currentLocation - 1][i]) * timep;
            });
        };

        String outputHeatFile    = dir + "heatFile.dat";
        String outputElasticFile = dir + "ElasticFile.dat";
        String outputStrainFile  = dir + "StringFile.dat";
        String outputStressFile  = dir + "StressFile.dat";
        String outputEqStressFile  = dir + "EqStressFile.dat";
        String deformFile =  dir + "DeformFile.dat";
        FES3T41 fsHeat = new FES3T41(mesh2C);
        FES3T42 fsElastic = new FES3T42(mesh2C);
        double kxx = 237, density = 2.7e3, specificHeat = 896.0;
        double elasticE = 6.9e10, elasticNu = 0.28e0, expansion = 2.36e-5;
        double beta = elasticE / (1 - 2 * elasticNu) * expansion;
        double theta = 0.5;
        SRMatrix LMatrix = new SRMatrix(fsHeat.GetNdof());
        SRMatrix RMatrix = new SRMatrix(fsHeat.getNdof());
        fsHeat.assembleHeatThetaScheme(new double[]{kxx, density * specificHeat, dt, theta}, new BVPType[]{BVPType.COMMON, BVPType.CONSIST_MASS}, LMatrix, RMatrix);
        DirectSSolver heatSolver = new DirectSSolver(LMatrix);
        heatSolver.Factorize_Cholesky(1);

        SRMatrix elasticK = new SRMatrix(fsElastic.GetNdof());
        fsElastic.assembleStiff(new double[]{elasticE, elasticNu}, BVPType.COMMON, elasticK);
        fsElastic.applyBC_MBN(elasticK, Direct.All, 10);
        DirectSSolver elasticSolver = new DirectSSolver(elasticK);
        elasticSolver.Factorize_Cholesky(1);

        double[] x = new double[fsHeat.getNdof()];
        double[] RHS = new double[fsHeat.getNdof()];
        double[] realTemp = x.clone();
        Arrays.setAll(realTemp,i -> x[i] + 300);
        setCurrentFlux.accept(time);
        Arrays.setAll(previousFlux, i->currentFlux[i]);
        double[] heat_loadFlux = new double[fsHeat.getNdof()];
        fsHeat.assembleFlux(new double[][]{currentFlux}, BVPType.COMMON, heat_loadFlux, 0);
        double[] heat_loadFlux0 = heat_loadFlux.clone();
        mesh2C.toTecplot(outputHeatFile, time, realTemp);
        double[] realHeatFlux = new double[fsHeat.getNdof()];
        double[][] vecX = new double[1][];

        double[] heatLoad = new double[fsElastic.GetNdof()];
        double[] totalLoad = new double[fsElastic.GetNdof()];
        double[] totalDisp = new double[fsElastic.GetNdof()];
        double[][] disp = new double[mesh2C.getnDim()][fsHeat.GetNdof()];
        double[][] strain = new double[6][fsHeat.getNdof()];
        double[][] stress = new double[6][fsHeat.getNdof()];
        double[] eqStress = new double[fsHeat.GetNdof()];
        mesh2C.toTecplot(outputElasticFile, time, disp);
        mesh2C.toTecplot(outputStrainFile, time, strain);
        mesh2C.toTecplot(outputStressFile, time, stress);
        mesh2C.toTecplot(outputStressFile, time, stress);
        mesh2C.toTecplot(outputEqStressFile, time, eqStress);
        while(time < endTime) {
            time += dt;
            setCurrentFlux.accept(time);
            System.out.println("Time: " + time);
            fsHeat.assembleFlux(new double[][]{currentFlux}, BVPType.COMMON, heat_loadFlux, 0);
            double[] Rx = RMatrix.mul(x);
            Arrays.setAll(RHS, i-> Rx[i] + 0.5 * dt * (heat_loadFlux[i] + heat_loadFlux0[i]));
            heatSolver.Solve_Cholesky(RHS, x, 1);
            Arrays.setAll(heat_loadFlux0, i -> heat_loadFlux[i]);
            Arrays.setAll(realTemp, i -> x[i] + 300);
            mesh2C.toTecplot(outputHeatFile, time, realTemp);

            setCurrentPressure.accept(time);
            double[] elastic_loadPressure = new double[fsElastic.GetNdof()];
            fsElastic.assemblePressure(new double[][]{currentPres}, BVPType.COMMON, elastic_loadPressure, 0);
            Arrays.fill(heatLoad, 0.0);
            vecX[0] = x;
            fsElastic.assembleHeatLoad(new double[]{beta}, vecX, BVPType.COMMON, heatLoad);
            Arrays.setAll(totalLoad, i-> heatLoad[i] + elastic_loadPressure[i]);
            fsElastic.applyBC_MBN(totalLoad, Direct.All, 0.0, 10);
            elasticSolver.Solve_Cholesky(totalLoad, totalDisp, 1);
            fsElastic.extractComponent(totalDisp, disp);

            mesh2C.toTecplot(outputElasticFile, time, disp);
            fsElastic.computeStrainStressEqStress(BVPType.COMMON, new double[]{elasticE, elasticNu}, totalDisp, strain, stress, eqStress);
            mesh2C.toTecplot(outputStrainFile, time, strain);
            mesh2C.toTecplot(outputStressFile, time, stress);
            mesh2C.toTecplot(outputEqStressFile, time, eqStress);
        }
    }
}
